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Abstract — The series expansion at the origin of the Airy 
function Ai(x) is alternating and hence problematic to evaluate for 
X > due to cancellation. Based on a method recently proposed 
by Gawronski, Miiller, and Reinhard, we exhibit two functions 
F and G, both with nonnegative Taylor expansions at the origin, 
such that Ai(x) = G(x)/F(x). The sums are now well-conditioned, 
but the Taylor coefficients of G turn out to obey an ill-conditioned 
three-term recurrence. We use the classical Miller algorithm to 
overcome this issue. We bound all errors and our implementation 
allows an arbitrary and certified accuracy, that can be used, e.g., 
for providing correct rounding in arbitrary precision. 

Index Terms — Special functions; algorithm; numerical evalu- 
ation; arbitrary precision; MiUer method; asymptotics; correct 
rounding; error bounds. 

Many mathematical functions (e.g., trigonometric functions, 
erf, Bessel functions) have a Taylor series of the form 

°° a" 
«=o 

For large jc > 0, the computation in finite precision arithmetic 
of such a sum is notoriously prone to catastrophic cancellation. 
Indeed, the terms \y„x"\ are first growing before the series 
"starts to converge" when n" > ax. In particular, when n'^ ^ ax, 
the terms y„x" usually get much larger than y{x). Eventually, 
their leading bits cancel out while lower-order bits that actually 
contribute to the first significant digits of the result get lost in 
roundoff errors. 

This cancellation phenomenon makes the direct computation 
by Taylor series impractical for large values of x. Often, the 
function y{x) admits an asymptotic expansion as x — » +oo that 
can be used very effectively to obtain numerical approximations 
when X is large, but might not provide enough accuracy (at 
least without resorting to sophisticated resummation methods) 
for intermediate values of x. 

In the case of the error function erf(x), a classical trick going 
back at least to Stegun and Zucker IfTTl is to compute erf(x) 
as G(x)/F{x) where F(x) - and |[l] Eq. 7.6.2] 

Gix) = erf(x) = ^ V , , (2) 
Vtt 1 . 3 . . . (2« + 1) 

The benefit of this transformation is that F and G are power 
series with nonnegative coefficients, and can thus be computed 
without cancellation. Algorithms based on Q tend to be 
competitive in some range a < x <b where x is large enough 
for cancellation to be problematic but small enough to make 
the use of asymptotic expansions at infinity inconvenient. Note 



that the obvious way to compute y{x) - e"* for x > fits into 
the same frame, now with G(x) - 1 and F{x) - e". 

Gawronski, Miiller and Reinhard [SJ, [14] provide elements 
to understand where these rewritings "come from". They relate 
the amount of cancellation in the summation of a series ([T]) to 
the shape of the Phragmen-Lindelof indicator of y, a classical 
tool from the theory of entire functions ifTOll . This description 
allows them to state criteria for choosing auxiliary series 
suitable for the evaluation of a given entire function in a 
given sector of the complex plane. They apply their method 
(called the "GMR method" in what follows) to obtain "reduced 
cancellation" evaluation algorithms for the error function and 
other related functions in various sectors. 

In this article, we are interested in the evaluation for 
positive X of the Airy function Ai LL Chap. 9]. The function 
Ai(x) satisfies the linear ordinary differential equation (LODE) 

Ai"(x)-xAi(x) = (3) 

with initial values 

Ai(0) = A := 3-^^^T{lr\ Ai'(O) = -B := -3"'/^''r(i)-'. 

The classical existence theorem for LODE with complex 
analytic coefficients implies that Ai(jc) is an entire function; 
and solving (|3]l by the method of power series yields the Taylor 
expansion Ai(x) = Af{x^) - Bxg(x^), where 

SI <3"'' ii 

Observe that while / and g are easy enough to evaluate 
individually, the difference Af{x^)-Bxg(x^) causes catastrophic 
cancellation when computed in approximate arithmetic. 

Using the GMR method, we derive a reduced cancellation 
algorithm for computing M{x). To our best knowledge, our 
algorithm for evaluating M(x) is new, and is the most efficient 
multiple-precision evaluation of Ai(jc) when x is neither too 
small nor too large, while the precision is not large enough to 
make methods based on binary splitting |4| competitive. 

Besides the new application, the main difference between the 
present article and the work of Gawronski et al. is our setting of 
multiple precision arithmetic "a la MPFR |7|". Specifically, on 
the one hand, we are interested in arbitrary precision arithmetic 
rather than machine precision only. This makes it impossible, 
for instance, to tabulate the coefficients of auxiliary functions 
when these turn out to be hard to compute. Also, we are 
looking for rigorous error bounds instead of experimental 
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Figure 1. The indicator function h of Ai. 

error estimates. On the other hand, we restrict ourselves to 
numerical evaluation on a half line instead of a complex sector. 

This article focuses on providing a complete algorithm in the 
specific case of Ai(x), x > 0. Yet, it should also be seen a case 
study, part of an effort to understand what the GMR method 
can bring in the context of multiple precision computation, 
and, perhaps more importantly, how general and systematic it 
can be made. We discuss this last point further in Sec. |VIII| 

The rest of the text is organized as follows. In Sec. |I] 
we use the GMR method to choose the auxiliary functions 
F and G. Then, in Sec. |II] and [III] we derive a few mathematical 
properties of these functions, including recurrences for their 
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series expansions and various bounds. Sections IV 
the details of our algorithm and its error analysis. Finally, 
in Sec. |VII| we briefly describe our implementation of the 
algorithm and present preliminary experimental results. 

I. The GMR method 

We now review the GMR method and apply it to obtain 
candidate auxiliary series for the evaluation of the Ai function. 
Since the method itself is not crucial for our results, we 
summarize it in intuitive terms and refer the reader to the 
original works [S], fl4j for more careful statements. 

The starting point of the GMR method is the following 
observation. Let y(z) = I]„>o ynZ" be an entire function. Assume 
that we have, in some intentionally vague sense. 



y(re"') * exp(/i(0)/^) 



(4) 



for large r. (To make things precise, we would assume that 
y has finite order p, and that h is its indicator function with 
respect to p 110] .) 

We consider the computation of y{z) in floating-point 
arithmetic using its series expansion. It is well-known IS) 
that, if the sum is performed in floating-point arithmetic of 
precision f, the relative error between y(z) and the computed 
sum is roughly given by (2„>o \ynz"\)/\y(z)\. The sum 2„>o bnz"! 
is larger than max„>o IjnZ"!, and usually of the same order of 
magnitude. Therefore, the number of significant digits "lost by 
cancellation" is roughly 

lg(max|y„z«|) - lg\y(z)\. 

Denote M{r) = sup|,|^^|y(z)| for all r > 0. Cauchy's formula 
implies max„(|y„|r") < M(r), and under a suitable version of 
hypothesis (|4]), one can actually show that max„(|y„|r") » M{r). 
Hence, the loss of precision by cancellation in the evaluation 
of y{re'^) is about 

M(r) 



Is 



\y(re'0)\ 



[imcixh)-h(6)]rP. 



For instance, when the all have the same complex argument, 
the maximum of h is reached for = 0, in accordance with 
the fact that the sum is optimally conditioned. 



4/3 ; 




/ 




2/3 

/ 






;r/3 2;r/3 ^ 


-2/3 





4/3' 




^'~^\/'^^' 




-2n/3 -n/3 


n/3 2n/3 " 


-2/3 





Figure 2. Indicator functions for F (left) and G (riglit). 

In the case of the Airy Ai function, the following asymptotic 
equivalent holds as z tends to complex infinity in any open 
sector that avoids the negative real axis [1 , Eq. 9.7.5]: 

exp(-fz3/2) 



Ai(z) ~ Ai(z) :-- 



(5) 



2 V^z'/4 

Additionally, Ai(x) is bounded for x < 0. Hence, we may take 
p = 3/2 and 

= -|cos(§6'), -7T<9<7r (6) 

(see Figure [T]l. The loss of precision is roughly proportional to 
1 + cos(|0). It is minimal in the directions of fastest growth 
9 - +|7r, and maximal for 9 = 0. 

If now two entire functions y and F both satisfy conditions 
of the form Q with the same p but different h (say hy and hf, 
respectively), we may expect that 

G(z) = F(z)y(z) * sxp([hy(9) + /if (0)]zO. (7) 
The GMR method consists in reducing the summation of the 
series y for z in some given sector to that of an auxiliary series 
F(z) and a modified series G(z) related by (|7]l. The value of y(z) 
is then recovered as G{z)/F{z). The auxiliary series is chosen, 
based on the shape of hy, so that both hf and /ig = hy + hp 
take values close to their maximum in the sector of interest. 

Gawronski et al. usually take F{z) = exp(az'^) and search for 
a value of a that makes (max he) -he as small as possible on a 
whole complex sector The choice of exponentials as auxiliary 
series is not appropriate in the case of Ai, since exp(z'') is an 
entire function of z only for integer p. 

However, as we are interested in one direction only, we can 
easily build a suitable auxiliary series from Ai itself. Indeed, 
we may "shift" the indicator function of Ai by In/ 3 to the left 
or to the right by changing z to j'^'z, where j = e^"'^^. (Note 
that this is not the same as changing 9 to 9 + in (|6|.) When 
we add such a shifted indicator to the original h^i, one of the 
humps of the curve cancels out with the valley in the middle. 

Using this idea, we set 

F(x) = Aiijx)Ai(f^x), G(jc) = F(x) Ai(x). (8) 

The indicator functions of F and G are pictured on Figure [2] 
Based on their shapes, we expect that both series are optimally 
conditioned on the positive real axis. We shall prove that this 
is indeed the case in the next two sections. 

This second method of constructing auxiliary series seems 
to be new, and applies to many cases. For instance, applying 
it to the error function leads to 

2x Y l-3---(2«-l) ^„ 2„ 



F(x) = -/ erf ix = 



(2n + l)! 



-2"x'", 
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G{x) = F{x) erf x = > 

n=0 



(4n + 2)! 



■(4«)4„^4„ 



a slightly worse alternative to (j2|. The advantage of ([2| comes 
from the fact that e*" is faster to evaluate than F(x). 

n. The auxiliary series F 

The functions F and G being chosen, we need to establish 
appropriate formulas to evaluate them, along with error bounds. 
Much of our analysis will be based on the following simple 
estimate fTT, Chap. 4, §4.1], where Ai is the leading term of 
the asymptotic expansion of Ai(z) at infinity from Eq. (|5]l. 

Lemma 1: The Airy function Ai satisfies 

|Ai(re''')/Ai(re'*) - 1| < j]i(0)r-^'^, \e\ < n, 

where 771(6!) = ^(cos jT^'^. 

Now consider the power series expansion of F at the origin: 



F(x)^Ai(jx)Ai(j-'x)^J]^"^"- 



(9) 



«=o 



Proposition 1: The coefficients F„ are positive and satisfy 
the two-term recurrence relation 

in + l)(n + 2){n + 3)F„+3 - ^iln + 1)F„ = (10) 

with initial values 

Fo = 3-^'^r(j)-^, Fi = (2 ^7t)-\ F2 = 3-^'^T(})-^. 

Proof: As a general fact, if two functions w and y each 
satisfy a linear homogeneous ODE with coefficients in Q(x), 
then their product wy satisfies an equation of the same class 
that can be explicitly computed IfTSl Sec. 6.4]. The functions 
Aiij'^^x) satisfy the same differential equation ([3]) as Ai(x) 
itself. (In other words, the substitution x y'^'jc leaves ([3]) 
unchanged.) Applying the procedure mentioned above to two 
copies of that equation yields F^^\x) - 4xF'{x) - 2F{x) = 0. 

Similarly, when an analytic function y satisfies a homoge- 
neous LODE over Q(x), we can compute a recurrence relation 
with coefficients in Q(«) on the coefficients y„ of its power 
series expansion. In the case of F, we get ( [TO] i. Finally, we 
compute the initial values F\,F2,Ft, from the first few terms 
of the Taylor expansion of Ai(x): 

F{x) ^(A-Bjx + 0(x^))(A - B f^x + 0(x^)) 

= - ^ABx + B^x^ + 0(x^). 

It is then apparent from ( fTO] ) that F„ > for all n e N. ■ 

Thus, the coefficients of F{x) obey a two-term recurrence 
whose coefficients do not vanish for n > 0. This allows one to 
compute them in a numerically stable way (see Sec. |VI| i. 

III. The modified series G 
Recall that G(x) - Ai(x) Ai(jx) Ai(j'^ x), and set 

exp(|x-'^^) 



G(x) = Ai(x)Ai(;x)Ai(7"'jc) 



(11) 
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Proposition 2: The function G is an entire function with 

power series expansion of the form G{z) - Yj7=Q^nZ'^" ■ The 

coefficient sequence (G„)„eN is determined from its first terms 

,1 , 1 1 

Go=a3 = — — ^, Gi = — -B^ 



9r(2/3)3' ' 2 18r(2/3)3 3r(l/3)3 

by the recurrence relation 

(n+l)(3n+4)(3n+5)(n+2)G„^2-l0(n+lfG„+i+G„^0. (12) 



Proof: First, observe that G(jz) = G{z) and G(z) - G(z), 
so that the Taylor expansion of G at the origin is a power 
series in with real coefficients. The same routine reasoning 
as in the proof of Prop. [T] yields the LODE 

G^*\x) - 10xG"(x) - lOG'(x) + 9x^G(x) = 0, 

and from there the recurrence ( [T2] i. The coefficients of ( [T2] l 
do not vanish for n > 0, so that the sequence G„ is indeed 
determined by Go, Gi, and ( [T2] l. ■ 
Nonzero solutions of ( [T2] l decrease roughly as n!"^ for 
large n. Setting c„ - nl^G„ yields the "normalized" recurrence 
(3n + 4)(3n + 5) 



-c„+2 - 10c„+i + c„ 



0. 



(13) 



(n + l)(n + 2) 

Letting n go to infinity in the coefficients of ( pj) , we get a 
limit recurrence with constant coefficients whose characteristic 
polynomial 9a^ - lOo- -1- 1 has two roots of distinct absolute 
value, namely ^ and 1. By the Perron-Kreuser theorem ||20l 
Theorem B.IO], it follows that any solution (v„)«eN of ( [T2| l 
satisfies v„+i/v„ ~ an^^ with either a = 1 or a = ^. Solutions 
(v„) such that v„+i/v„ ~ ^ are called minimal and form a 
subvector space of dimension 1 of the solutions of ( [T2] i. 

It turns out that (G„) is a minimal solution of ( [T2] i. But our 
analysis uses a bit more than the rough estimate G„ ~ nl^^9^". 
Prop. [3] below provides a more precise estimate which implies 
the minimality. Before turning to it, we recall a standard bound 
on the tails of incomplete Gaussian integrals |1 Eq. 7.12.1] 
and state a second technical lemma. 

Lemma 2: The complementary error function erfc x = 1 - 



erf X = ^ r°° e ''df satisfies < erfcx < 

^J7: Jx 

Lemma 3: The expression 

-;r/3 



e ^' for X > 0. 



/(r) 



f 



;,J'2(-l+COsf) 



satisfies < I(r) < 0.51 for all r > 10. 



Proof: We first use the inequality cos(|6') < 1 ^ 



•■2"' - ^ iO" 
(valid for < 6 < n/3) followed by the change of variable 

ip = ^|3f5r'^l^ e io get 

I(r) < r'/' e-i^'^'^de = ^ r''' e'^'d^, 

Jr-^l« V3 Ja{r) 

where a(r) = 



V375ri^^ and b(r) = 
r^'d^o < I e-'^'d(fi 

Ja(r) 



-^r^^'^. Now we have 



2 

_3 1/4 



erfc 



(^/3r'' 



V5 

and, by Lemma [2] I(r) < |r'^'*e s*^"^. This last bound is 
decreasing and less than 0.51 when r = 10. ■ 

The following proposition is the main result of this section 
and the starting point of much of the error analysis that follows. 
Though somewhat technical, the proof is mostly routine. 

Proposition 3: The sequence (G„) satisfies 
^ 1 

G„~r„ = — -, « ^ 00, 

4 V37r9"n!2 

with relative approximation error \G„/y„ - 1| < 2.4n"'^"*. 

Notation 1: We write ^(ei, . . . ,e„) - 20^/c|i,n] Wiei^i^ so 
that, if \0i\ < Si for all then nf=i(l + 0,) = 1 + with 
\0\ < E{e\, . . . ,e„). Note that we also have E(asi, . . . ,as„) < 
a E{ei ,...,£„) for < or < 1 . 



Proof: The proof is a standard application of the saddle 
point method |6, §VIII.3], which we work out in some detail 
in order to get an explicit error bound. 

Start with the estimate from Lemma [I] Writing z - re'" with 
r > 1 and \6\ < n/3, it follows that 

\G(z)/G(z)-l\<m(d)r-''^ (14) 
where 772(6!) = £(771(6), 771 (6* + f), 771(61- f )). Note that 772(61) 
increases with \0\. Using the connection formula |[1] Eq. 9.2.12] 

Ai(z) + jAi(jz) + j-'Ai(j-'z)^0, 

and assuming now that < < 7r/3, we see that \G(z)/G(z)\ is 
bounded by 



1 Ai(z) 1 


|Ai(r'z)| 


1 Aife) 1 


|Ai(z) 


1 |AiO-'z)| 


|AiO-'z)| 


1 Ai(z) 1 


lAiO-'<-)l 


U AiUz) 1 


lAi(z) 


1 AiOz) 1 


\Ai(j-h)\ 



<(l+77i(0)r-^'O(l+ 771(0- f)r-3/2) 

(e-3™^¥(l +77i(6l)r-5)+ 1 +771(0- f )r-2). 
It follows that 

\Giz)\ < 2A\G(z)\, \z\ > 10. (15) 
By symmetry, these inequalities hold for -7r/3 < 6 <0 too. 
Now fix n > 10 and write 



1 

27r7 



G(z) 

,3n+l 



dz: 



3 G(r) , 



p/3 



G(re'") _ 



Z-'-' 2n r^" ' ' J_^/3 G{r) 
As prescribed by the saddle point method, we choose 

r = (3« + 3/4f'\ 

so that 

G(z) _ G(r) o((,_,)2) 



(i.e., £ log(G(z)z-3")|,=. = 0), and 

-5/8 



-5/12 



yo = r = (3n + 3/4)" 

Observe that r > 10 and < 6/0 < 1/4 < 7r/3. 
With this choice of r, we have 

V3 G(r) n\^e^" Ve 



(16) 
(17) 

(18) 



Stirling's formula in the form IfTSl 
combined with the bound 



n^"27m (I + ^) 



I n" V2;r« 



L^2n+1 ■ 

- 1| < s(n) = e li; - 1, 



(MjmF - li' yields 



I// - II <£(.(«), .(«),^,i)<l, «>11. 



(19) 



We further let 



h 
h 



G(re'") _ 
G{r) 



J-e, 



3m0 



2V^ _ 



3/4 



r'/^0^\d6. 



and now aim at bounding the relative difference I/1//4 - 1| 
First, using ([TSj, we have 



G(r) 



-d9 < 4.2 



Gir) 



-do 



-;r/3 



= 4.2 



r\i'-"^(-'--¥)d0 



and hence 
by Lemma [3] 



|/i-/2l<2.15r 



-9/8 



(20) 



Next, thanks to ( [T7] i, we can write 

Gire-'ye-'"'" = (5(r)ei^"'(«''P(3'fl)-i-i«) = G(r)e(-i»'+"(''»^"' 
where 

\u{e)\ = f |e5* - 1 - lie + |6»2| < l9^. (21) 

Let 

G(re''')e-3'"® 



viO) 



G(r)exp(-|r3/26i2) 



1. 



Since fl^r^^^ - r by the choice ( [T8] l, using ( |2T] i and the 
inequality |e~ - 1| < |z|e^'', we have 

|/'Vfl)_i|<!,-3/8^^-'^<0.44r-3/8, 
8 

for |6i| < 00 and r > 10. By ([T4]i, it follows that 
G(r£^ 3,2„(g) 



|y(0)| = 



- 1 



G(re'») 

< £(772(0o)'-"^^^O.44r-^'"^) < 0.95r-^'\ 



We get 



1/2-/31 = 



J-e, 



e-'^'"''>\i0)d0 



< 0.95r-^^^h 



(22) 



Finally, is an incomplete Gaussian integral: 

y3^3/4 

Lemma [2] yields 

2r-i/8e-J'-"' ,,„ . 

I/3//4 - 1| < ^ < 0.66r-'/^e- 



37T 



(23) 



Putting together ([19] |20 

|/i 



I/1//4-II < 



■/2I 



22] |23l), we obtain 
1 



h hh 

< LOer"^/*^ +£(0.95r"^''^0.66r" 

< 2.45r-3/8 < L9n-'/^ 



It follows that 
G„ 

7n 



G„ ^ 



h 



< E(0.5n-\ 1.9«"'^'*) < 2.4«"'/'* 



for 71 > 11. It is easy to check that the final bound is valid for 
1 < 77 < 10 as well. ■ 

Note that the above bound is not the best we can get by this 
method. Indeed, by choosing the exponent of r in ([TSj closer 
to -3/4, we obtain |G„/7„ - 1| = 0{n'^^^ + e) for any e > 0. 
This comes at the price of a larger constant factor and thus 
more terms to check separately. 

A first consequence of Prop. [3] is that the series G{x) has 
nonnegative coeflicients, as stated below. We also deduce 
several other technical results that will be used to bound various 
error terms. Recall that c„ - n\^G„, and let r = 3/20. 

Corollary 1: The sequence (c„) satisfies < c„+i/c„ < r for 
aU 72 > 0. Accordingly, we have < G„+i/G„ < r(77 + 1)-^. In 
particular, the G„ are positive. 

Proof: Prop. [3] implies that c„ = (1 + 0„)/(4 V3;t9") with 
Wnl < 2An-^l\ When tj > tj,, = 70000, we have |0„| < 0.148. 
This implies c„ > and c„+i/c„ < (1/9)(1 + 0„+i)/(l + 0„) < r. 
The inequalities up to tj = tjq are checked by computing the 
corresponding terms with interval arithmetic. ■ 

Corollary 2: For any 7? > 1, it holds that G„ < (e/(377))^". 



Proof: With 9„ as in Corollary [T| we have G„ < y„(l +0„), 
and hence G„ < l/(9"n!2) < (l/9)"(e/n)2«. ■ 
Lemma 4: When x > and + 1 > y/lrx^^^, we have 

N N-l 

Proof: The first inequality is obvious as G„ > 0. To show 
the second one, observe that a fortiori Tx^/{n + 1)^ < j for 
all n> N. Using Corollary 111 we deduce G„+\x^/G„ < j, and 
hence Zn>N G„x^" < Gnx^^I + i + i + ■ ■ ■ ) = 2Gnx^^. ■ 
Lemma 5: The function G satisfies 



Algorithm 1: Evaluation of G 



for all jc > 1/2. 

Proof: Let n(x) - E(rii(0)x^^^^, crx'^^^, crx'^^^) where cr - 
/7i(f ) = I V2. Lemma [l] implies |G(x)/G(x) - 1| < ju(x) for 
X > 0. We can check that fi(x) is a decreasing function and 



0.01 < 



1-M3) 



1 +fi{3) 



< 0.04, 



87r3/2 ' 8;r3/2 
whence the desired inequality for x > 3. For 0.5 < x < 3, using 
the bounds < e'' - (1 + x + jx^ + ^x^ + j^x'*) < ^x'* (valid 
for < X < 1 3^^^) and Lemma |4j we are reduced to checking 
explicit polynomial inequalities. ■ 

IV. Roundoff error analysis 

We now turn to the floating-point implementation of the 
functions F{x) and G(x). To make the algorithm rigorous, we 
will use classical techniques of error analysis that we briefly 
recall here. We refer the reader, e.g., to Higham ||9l, for proofs 
and complement of information. 

Notation 2: For x 0, we denote Exp(x) = Llog2 \x\i + 1, so 
that 2^"?^-' <\x\< 2'^''pW. 

Notation 3: If x e R, o(x) denotes the floating-point number 
closest to X (ties can be decided either way). Circled operators 
such as © denote correctly rounded floating-point operations. 
Floating-point numbers are represented with t bits of precision. 

Thus, we have o(x) = x{l + 6) and x - o(x)(l + 6') with 
< 2"'. We will also extensively use the relative error 
counter notation z (k). 

Notation 4: We write T = z{k} when there exist 6 1, ... ,6k 
such thatT= znti(l + ^O'^' with < 2"' for all /. 

Roughly speaking, each arithmetical operation adds one 
to the relative error counter of a variable. The overall error 
corresponding to an error counter can be bounded as follows. 

Proposition 4: Suppose that we can write z{k) and that 
k2-' < 1/2. Then?= z(l + 0) with \6\ < 2k ■ 2"'. 

V. Evaluation of the modified series 

As we shall see in the next section, evaluating the auxiliary 
function F is fairly straightforward. The evaluation of G is 
more involved. Indeed, while 2 G„x^" is well-conditioned as a 
sum for X > (this is the whole point of the GMR method), the 
minimality of the sequence (G„) among the solutions of ([12]) 
implies that its direct recursive computation from the initial 
values Go and Gi is numerically unstable (cf. EOl ). 

There is a standard tool to handle this situation, namely 
Miller's backward recurrence method O, l20ll . Miller's method 



Input: a target precision > 1, a point x > 0.5 
Output: s such that |G(x) - i| < 3 ■ 2-p G(x) 



-lr-3/2 



1 Set a,j3, 6, y s.t. a < 3 e 
r > l/log2(20/3), 8 > (2/3) log2(e)((20/3)'/2 _ l)x 



/?<(2/3)log2(e)x 

1/2 . 



3/2 
//2 



2 N) ^max(l,r(3/10)'/2x 



1/2 3/2 . 



11); 



Exp(x) - [j6J; 



3 Set > A^o s.t. Exp((aA^)2'^) >p + 9+^ 

4 SetR> maxiN, yip + 2 + 5)); 

5 Set t s.t. 128 (A^ + 3) 2"' < 2''' and (R + 2) 2'' < 2"^; 

6 for (fl ^ 1, <- 0, / ^ - 1; / > 0; / ^ / - 1) do 

7 c <— a; 

8 fl ^ (10 (8> fl) e (3/ + 4)(3i + 5)®b0 ((/ + 1)(/ + 2)); 

9 b <— c; 

10 if / = A^ - 1 then s <— a; 

11 else if / < A^ - 1 then 5 <— a © (x-' ® 5 (/ -i- 1)^); 

12 return (s0a)0 o(9 r(2/3)^); 

allows one to accurately evaluate the minimal solution (c„) of 
a recurrence of the form 

fl2(n) Mn+2 + fli(n) Mh+1 + flo(n) «« — 0, ao(n)a2(n) + 0. (24) 
The idea is as follows: choose a starting index R and let 
(arbitrarily) - and = 1. Then compute m„ as 

-flo(n) ' (fli(n) Mn+i + fl2(«) M„+2), n — R — 1,...,1,0. 
It turns out that, for large R, the computed sequence (m„) is 
close to a minimal solution of the forward recurrence. Since all 
minimal solutions are proportional to each other, we recover 
an approximation of c„ as c„ » {co/uo)u„. 

We use Miller's method to evaluate the minimal solution (c„) 
of the normalized recurrence ( [T3] l, and we get an approximation 
of G(x) = E«>oG„x3" ^ (co/mo) I^^^o u„x^"/in\^). The 

algorithm is summed up as Algorithm [1] The rest of this section 
is devoted to its proof of correctness, i.e., the proof that the 
value s it returns satisfies |G(x) - i| < 3 ■ 2^pG(x). 

Proposition 5: With A^ as in Algorithm [1] the truncation 
error satisfies | ZZnGhx'^"] < 2-p G{x) for all x > 0.5. 

Proof: First, because of Line [2] of the algorithm, we have 



■^x^ < (A^-H 1)^, so that 



W.G,. 



n=N 



< 2Gnx^^ < 2 



jj(.3/2 



2N 



by Lemma |4] and Corollary |2] Line |3] then ensures that 

r3/2N2A' 



< 2 iaN) 



-IN 



2-P 



-3/4 



< 2-''G(x), 



3N j ~ ~ 128 

the last inequality coming from Lemma [5] ■ 
There are two sources of error besides the truncation: first, 
(m„) is not exactly proportional to (c„), especially when n is 
close to R. Second, roundofi" errors happen during the evaluation 
of (m„). Rigorous bounds for both sources of error have been 
proposed by Mattheij and van der Sluis |11|. We combine 
them with classical techniques (well-explained, e.g., in |5|) to 
choose the starting index R and the working precision t so as 
to to guarantee the final accuracy. 

We now recall Mattheij and van der Sluis' main result 
(adapted to our particular case, which simplifies the statement 



quite a bit). Consider a recurrence of the form (|24]). Denote 
by (c„) a minimal solution that we wish to evaluate, and let 
(d„) be the solution such that do = di - 1. Assume that (d,,) 
is a dominant solution and that the sequences (c„), (d„) and 
(c„/d„) are decreasing. Define H - ^ Tjh^o j. ^nd, for / < R, 

= | ■ ^ i,^, landB, = 



-ai(') 
ao(') 
1 



-aiji) 
ao(0 





Let vj? e be a column vector, and for /< 7? - 1, let v; be the 
result of the floating-point evaluation of B, v',+i at precision t. 
Write v; = (m,, m,+i)^. If all operations were exact, (m,) would 
be the solution of the recurrence such that {ur, mr+i)^ = vr. To 
take rounding errors into account, we write v, = (B, +2"'g',)v,+i 



Altogether, this ensures that < 1. So, writing = 
1 + 6i, we have 

< 2\^JLi\ + |6l,| < 1 U.5{N + 3)2"' + T^-' < 2-P + t*"', 

and therefore \s - ^^q' GiX^'\ < 2-pG(x) + t'^G(t-^'^x). 
Lemma 6: For x > 0.5, we have G{t^^^^x) < 2^''G(x). 
Proof: It follows from Lemma [5] (and r < 1) that 

> 1 expf ^.3/2(1 _ ^-1/2) 

G(t-i/3x) 4 ^ 
and the algorithm ensures that 



R > 



1 



logjT 



p + 2+-log2(e)(T-'/2. 



3/2 



for some matrix ^, instead. Define = (3'Ri,yR2)^ = ^^r'^R- whence < i 2 ^' exp(| .^3/2(1 _ r i/2)) ^nd the result 



Let "Fi - \\U^ 0iUi+\\\, the matrix norm being subordinate to 
the II ■ lloo norm for vectors, and let ;F > max, ;F/. 

Theorem 1: IfTTl Theorem 4.1] Provided that the quantities 

CRdoyui \\yR\U 



TR2- 



■(R + H){\3T2-') 



dRCoym lyml 
are all bounded by 0.1, the approximate value m,- computed by 
Miller's algorithm satisfies (co/mo)m, - c, (l + O,) for some 0, 
such that \6i\ < T, + Ri, where 



Ti = 1.5 



cr di cr do 



codR 



JRI 



yR\ 



Ri = 1.5 



yR\ 



s(i + 2H). 



[cidR 

and e= 1.3:r2-'. 

Turning back to the special case c„ = n\^G„, Theorem [T| 
applied to ([T3| yields the following. Recall that t - 3/20. 

Corollary 3: Set vr = (1,0)^ and 



10 
1 



-r(i) 




where r(n) - 



(3n + 4)(3n + 5) 



(n + l)(n + 2) 

Then, in the notation of Theorem [l] we have T,- < ' and 
Ri < 76.5(/ + 4)2-' < 76.5(N + 3)2"' for all / < N. 

Proving this corollary still requires some work. We postpone 
it for a bit to explore the consequences of this statement. 

Observe that lines |8] and |9] of Algorithm [T| are equivalent 
to a floating-point evaluation of B,v,+i where v,+i - {a,b)^. 
Hence, at each loop turn, we have a - u, just after line |8] 
Lines [10] and [TT] are a Horner-like evaluation scheme, so that 
s * Z/Io' UjX^'/i]^ at the end of the loop. More precisely, 
assuming x^ is approximated by o(x) ® o(x) ® o(x) and division 
by (; + 1)^ is performed as two successive divisions by (/ -i- 1), 
an easy induction shows that one can write 

A'-l 

s = yiu,^\{9(i+i)}. 



JV — i / 

/■=() ^ 



,-l2 



Line[T2]adds 3 to all error counters. Using Prop.[4]we conclude 
that, at the end of the algorithm, we have 



A'-l 



A'-l 



where |/i,| < 2 (9(; + 1)-h3)-2-' < 18(A^ -1-3)2-' and |0,| < Ti+R,. 

Since R > N, we have T, < t for all / < A^. The definition 
of t implies 128(A^ -i- 3)2"' < 2-'' < 1/2, hence < 76.5/256. 



Since \G(x) - s\ < \G(x) - Y^to G/x^'l + \s- ^f^o' Gix''\, the 
correctness of Algorithm [T] is a consequence of Prop. [5] and 
the above discussion. 

Theorem 2: The value s returned by Algorithm [1] satisfies 
\G(x) -s\<3- 2-PG(x). 

The remainder of this section is devoted to the proof of 
Corollary |3] We begin with a crucial lemma. Let {d„) be the 
solution of ( [TJI defined hy do = di = I. 

Lemma 7: For all « > 1, we have d„^\ < d„ < (1 +ri{n))d„+i 
with T](n) = l/(3«2). 

Proof: We proceed by induction. Since d2 = 9/10, the 
property is true for « = 1. Now, supposing it for an arbitrary «, 
we get (9 - T]{n))d„+i < lQd„+i - d,, < 9d„+i, so 



9 - i](n) 



dn+l ^ dn 



+2 



9 ^ 
r(n) 



rin) = 



(3n + 4)(3n + 5) 



r{n) ~ ~ r(n) (n + l)(n + 2) 

We conclude by observing that 9 / r{ri) < 1 and r(n) /(9 - ?/(«)) < 
I +ri(n + 1) for n > 1. ■ 
Corollary 4: For all «, we have 0.783 < d„ < I. 
Proof: By Lemma [7] (d„) is decreasing and do - 1: this 
proves the right-hand side. For n > 100, Lemma |7] implies 

n-l it 

-^100 < n ''•^'^^ - " n*^^ 

1=100 P^'^ 1=1 

Using the exact value poa - ^ sinh(^) [1, Eq. 4.36.1], we 
check that d„ > dioo pgg > 0.783 for n > 100. As {d„) is 
decreasing, the inequality holds for n < 100 too. ■ 

This estimate, combined with Corollary [T] gives almost all 
we need to check the hypotheses of Theorem [T[ We use the 
notation introduced for the statement of Theorem [T] 

Corollary 5: The sequences (c„), {d„) and {c„/d„) are de- 
creasing. Moreover, the following inequalities hold: H < 2, 



\yR2ly, 



fill 



IbslU/b. 



R\\ 



1, and dett/, > |c^. 



Proof: Corollary [T| shows that (c„) is decreasing and 
Lemma |7] shows that {d„) is decreasing. Together, they imply 
TtJt - + l/3n^) < 5 for n > 1, and we can check that 
c\ld\ < Col do- Hence (c„/d„) is decreasing. 

Corollary [T] also implies c„ < t"co for any n, and Corollary |4] 
shows that do < (i„/0.783 for any n. It follows that 

«-i 

- > ^ < :r^^T' <2. 



0.783 



Wehavec/oM = l,dyld2 = 10/9, and c/^M+i < 1 + 1/(3/^) < 
10/9 for / > 2, and, by definition of Ur and v^, 

cr I dR+i/dR -1\/1\ cr IdR+l/dR 
-cr+i/cr 1 /\0/ dett/ft \-Cfi+i/cR 

< g. This implies ly^al < \yR] 



yR 



Ur\'R 



detC/R 

It follows that \yR2lyRi\ < f i 
and therefore HyfilU/lyRil = 1. 

Finally, from the expression det f/,- - cj{di+i/di - c,+i/c,), 
we obtain c7^ det t// > ( - r) = | . ■ 
Lemma 8: A suitable value for 'F is 39. 

Proof: The multiplication B, v, is performed on lines [8] 
and |9] Denoting by a' and fe' the new values of a and b, we 
can write b' - a and, with r(/) as in Corollary [3] 

a' = 10a <2) - r(i> (5) = a' = 10fl(l + 0) - r(iXl + 6') 

where |0| < 4-2"' and l^'l < 10-2"' by Prop.|4] (This assumes that 
the multiplication by r{i) is performed through four successive 
multiplications and divisions). Therefore, we have 

where, for matrices, \A\ < \B\ means that the inequality holds 
entrywise. Since moreover 

■c^ld^ldi 1^ 
{cMici lj-3c,i3/20 1, 

c,+2/c,+i di^2ldi^}] - (3/20 
and5^-<|f/ri|-|^,Mt/;+il, we get 



1^' l^det^/,.1 

|t//+ll < C;+l 



3 

20 



40+ |r(0 
l)(40+iK0) 



40 + 10r(0 
I5 (40 + lOr(O) 



Hence \\Ti\\ ^ 8 + ^ r(/) + 8 + 2 r(0). The result follows because 
r(i) < 10 for all / > 0. ■ 

We can finally prove Corollary [3] 

Proof: In order to apply Theorem [T[ we need to check that 
(i) 397;2-' < 0.1, (ii) ^ < 0.1 g^, and (iii) (R+2)(5Q.l ■!-') < 
0.1. By definition of t, R + 2 < 2'"^, so (iii) is satisfied. We 
notice that (iii) trivially implies (i). Finally, by Corollary |5] it 
suffices to prove < 0.6 to prove (ii). This is clearly true 
since doldR < 1/0.783 and cr/cq < < < t. 

In conclusion, we can apply Theorem [1] hence 



Ti = 1.5 



Cr di Cr do 



CidR 



codR 



yR2 



yRi 



1.5 

0.783 



Therefore, T, < 0.639 ■ r'' ' < t* ' as announced. Theorem [T| 
together with Corollary [s] give = 1.5 (50.7 ■ 2 ')(' + 2H) < 
76.5 ■ (/ + 4)2'. ■ 

VI. Evaluation of the auxmary series 

The implementation of the auxiliary series F is much easier 
than that of G: the recurrence is a simple two- term recurrence 
that can be accurately evaluated by forward recursion. We limit 
ourselves to a sketch of the (fairly standard) algorithm. 

A variable ao is used to successively evaluate Fq, F^x^, F^x^, 
etc., using the recurrence Q. Accordingly, two variables oi 
and 02 are used to evaluate the successive values of F-i^^ix^"^"^ 
and FT,n+2X^"^^- Each step adds 10 to the relative error counter 
of each variable. A variable s is used to accumulate the sum 



as the variables are updated. Therefore, after step K, we 
can write s = ^m)"' F,Jt:'<l + 10(//3) + 3K), or, bounding all 
errors uniformly 

N-\ 

i = ^F,x'<5(A^+ 1)), N^3K. 

1=0 

It is easy to see that q{ri) - 2(2n + l)/((n + l)(n + 2)(n + 3)) 
decreases for n > 1, hence the loop can be stopped as soon 
as (i) q(n)x^ < 1/2, and (ii) flo,fli,fl2 < 2^''p('>-''-1 These 
conditions ensure that the remainder of the series is bounded 
by (ao +ai+ 02) < ^2-''2E='pW. 

It is clear from Prop. Q that F„+3 < 4F„/n^ for any n, 
hence (using n! x (n/e)") we have F„ ^ (4e^/n^)"^-'. This is 
most likely an overestimation of the true value. Moreover, we 
can approximate F(x) by ^x^^^^ expi^x^^^) for x > 0.5 and 
by 2"^ + 2^*x for < jc < 0.5. These estimates are used to 
get a rough overestimation of the necessary truncation rank. 
Based on this estimation the working precision t is chosen so 
that 10(N + 1)2'-' < 2-^-P. Hence, we have | Zfjo ^ix' - s\ < 
2^^^P s < 2^''P(-')~^"'' . 

The initial estimation of the truncation rank is very unlikely 
to be underestimated. In the case it would be smaller than the 
actual truncation rank decided on-the-fly by the above criterion, 
this is checked a posteriori and, if necessary, the evaluation is 
run again with an updated working precision. 

We conclude by writing \F(x) - s\ < \F(x) - YjfJo^ Fi^'l + k - 
XfJo Fix'l < A2-P2E''P('') < 2-Ps. 

VII. Complete algorithm 

Assume p > 3. From the previous sections, it appears that we 
computed G such that G = G(x)(l+ei) with 16*1 1 < 3-2^'', and we 
coniputed F such that F(x) - F(l + 62) with |6'2| < 2"''. Hence 
GIF = {G{x)IF{x)){\ + 6i) where I6I3I = |6ii +02(1 +6li)| < 5-2-^. 
Indeed, the division is performed in floating-point arithmetic 
at precision p, leading a final result A — {G/F){1 + 64) with 
|04l < 2-P. Hence, finally, A = Ai(jc)(l + 65) with |0g| = |03 + 
6*4(1 +6*3)1 <7-2-'' <2-^J'-^K 

We developed a prototype implementation of this algorithm, 
based on the multiple precision floating-point library MPFR |7|. 
In order to test the implementation, we compared the results 
with that of the built-in implementation of Ai(x) available 
in MPFR, run with a larger precision. Random tests using 
several thousands of points x and target precisions p yield 
relative errors smaller than 2 between the result of our 
implementation and the result of MPFR, as predicted by theory. 

We also ran our implementation and MPFR with the 
same accuracy, in order to compare the performance of both 
implementations. The result is shown on Figure [3] (MPFR 3.1.1, 
GMP 4.3.1, on an Intel Xeon at 2.67GHz). Our method is 
faster than MPFR when x is fairly large. This is all the benefit 
of reducing the cancellation: when x is large, the working 
precision of MPFR has a large overhead because of the bad 
condition number of the series. In contrast, our method keeps 
a reasonable precision, of the same order of magnitude as 
the target precision. Both implementations are comparable for 
small values x and target precisions; and MPFR wins over our 
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Figure 3. Fastest method as a function of x and p (scales are logaritlimic). 

implementation only for intermediate values of x and fairly 
large precisions, when our implementation pays the cost of 
evaluating two series instead of a single one, with little benefit 
in terms of working precision. 

VIII. Outlook: Towards a GMR algorithm! 

As mentioned in the introduction, we see the present work as 
a case study. Indeed, in spite of the many technical details that 
occupy much of the space of this article, we really used very 
few specific properties of the Airy function besides Equation (|3]l. 
Looking back, our analysis essentially relies on the following 
ingredients. 

(i) The ability to find auxiliary series. The indicator functions 
used in the GMR method depend only on the behaviour at 
(complex) infinity of the entire functions they are associated to. 
In the case of the solution of a LODE with analytic coefficients, 
this behaviour is entirely determined by the differential equation 
along with a finite number of "asymptotic initial values" |19|, 
ifTSll . Once the indicator function is known, it remains to exhibit 
appropriate auxiliary series. Doing this in a truly general way 
remains an open problem. Yet, both the original GMR method 
and our variant apply to many cases, and it is likely that they 
can be combined and further generalized. 

(ii) An efficient way to compute their coefficients. This 
seems to be considered a major limitation in the original 
GMR paper |8|. But, as already mentioned, recurrences with 
polynomial coefficients automatically exist as soon as both the 
original function to evaluate and the auxiliary series satisfy 
differential equations with polynomial coefficients. Numerical 
stability is not much of an issue in the case of three-term 
recurrences, thanks to Miller's method, though many technical 
details must be settled. The situation is more complicated in 
general for recurrences of higher order Observe, though, that 
we proved the minimality of (G„) by using essentially the same 
asymptotic properties that were exploited by the GMR method 
in the first place. The minimality may hence not be fortuitous 
and might generalize. 

(iii) Upper and lower bounds on the coefficients and sums 
of the series F and G. All these bounds were derived, in a 
pretty systematic way, from the asymptotic expansion of Ai at 
infinity completed with Lemma [T] Bounds similar to that from 
Lemma [1] can themselves often be obtained from a LODE 1 13 1. 

(iv) Roundoff error analyses. Our error analyses follow a 
very regular pattern and could probably be abstracted to a more 



general case or automated. 

In short, most steps of our study could apparently be 
performed in a systematic, if not algorithmic, way starting 
from Eq. ([3]l plus a moderate amount of additional information. 
Systematizing the GMR method based on this observation 
seems a promising line of research. Solutions of LODE with 
polynomial coefficients are known in Computer Algebra as D- 
finite, or holonomic, functions. Powerful symbolic algorithms 
are available to work with this class of functions, and already 
cover points (i)-(iii) above to some extent Q. It is perhaps 
too ambitious, though, to aim to handle the whole class of 
D-finite entire functions. Thus, at this point, we suggest looking 
for a natural subclass to which the method applies in a truly 
systematic way, hopefully to the point that it can be automated 
using Computer Algebra. 

Acknowledgments: We thank Paul Zimmermann for point- 
ing out the GMR article to us. 
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